Lab 10: Multiple Regression

library(tidyverse)
library(openintro)
download.file("http://www.openintro.org/stat/data/evals.RData", destfile = "evals.RData")
load("evals.RData")

Exercise 1

This is an observational study in which the student responses are being observed. Since this is not an expiriment, there is no control group to compare the results too. So a better question would be:

“Does beauty have an impact on course evauluations, or is any impact just due to chance?”

Exercise 2

The score variable is centered around 4.6 with a right-skew and average variability and spread. The skew tells us that most students give above average ratings. I expected the distribution to be a little more centered around average ratings, but it makes sense that most students would rate favorably if they are rating at all.

hist(evals$score)

Exercise 3

The scatter plot of age and average beauty rating shows a negative linear trend with a moderately strong correlation. As age increases, the average beauty rating decreases linearly.

Explanatory Variable: Age Response Variable: Average beauty rating

plot(evals$age, evals$bty_avg)

Exercise 4

Some of the observations overlapped, but with the noise it is easier to see the true number of observations.

plot(evals$score ~ evals$bty_avg)

plot(jitter(evals$score) ~ evals$bty_avg)

Exercise 5

The equation for the regression line is: y_hat = 3.88034 + 0.06664*bty_avg

The score increases by 0.06664 points as the average beauty increases by a point. With an R-squared value of just 0.03502, there is a 3.5% chance of the results not being due to chance.

The correlation is also not very strong at 0.187.

Overall, it does not appear to be a practically significant predictor.

m_bty <- lm(score ~ bty_avg, data = evals)

plot(jitter(evals$score) ~ evals$bty_avg)
abline(m_bty)

summary(m_bty)

## 
## Call:
## lm(formula = score ~ bty_avg, data = evals)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.9246 -0.3690  0.1420  0.3977  0.9309 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.88034    0.07614   50.96  < 2e-16 ***
## bty_avg      0.06664    0.01629    4.09 5.08e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5348 on 461 degrees of freedom
## Multiple R-squared:  0.03502,    Adjusted R-squared:  0.03293 
## F-statistic: 16.73 on 1 and 461 DF,  p-value: 5.083e-05

cor(evals$bty_avg, evals$score)

## [1] 0.1871424

Exercise 6

Linearity: The data barely shows a linear trend with a weak correlation coefficient.

Nearly Normal Residuals: The histogram of the residuals resembles a normal distribution closely, with a slight right skew.

Constant Variability: The residual plot shows constant variability. The points on the normal probability plot vary a more near the larger x values than for smaller ones.

Independent Observations: We can assume from the data set that each observation is independent.

plot(m_bty$residuals ~ evals$bty_avg)
abline(h = 0, lty = 3)  # adds a horizontal dashed line at y = 0

hist(m_bty$residuals)

qqnorm(m_bty$residuals)
qqline(m_bty$residuals)  # adds diagonal line to the normal prob plot

Exercise 7

plot(evals$bty_avg ~ evals$bty_f1lower)

cor(evals$bty_avg, evals$bty_f1lower)

## [1] 0.8439112

plot(evals[,13:19])

m_bty_gen <- lm(score ~ bty_avg + gender, data = evals)
summary(m_bty_gen)

## 
## Call:
## lm(formula = score ~ bty_avg + gender, data = evals)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.8305 -0.3625  0.1055  0.4213  0.9314 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.74734    0.08466  44.266  < 2e-16 ***
## bty_avg      0.07416    0.01625   4.563 6.48e-06 ***
## gendermale   0.17239    0.05022   3.433 0.000652 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5287 on 460 degrees of freedom
## Multiple R-squared:  0.05912,    Adjusted R-squared:  0.05503 
## F-statistic: 14.45 on 2 and 460 DF,  p-value: 8.177e-07

Check for outliers: The histogram shows few to non.

hist(m_bty_gen$residuals)

There are no obvious patterns against fitted values.

plot(abs(m_bty_gen$residuals) ~ evals$bty_f1lower)

There are no obvious patterns when viewing the data in order.

plot(m_bty_gen$residuals)
abline(h = 0, lty = 3)  # adds a horizontal dashed line at y = 0

The variability does not change much between ethnicity and gender.

plot(m_bty_gen$residuals ~ evals$ethnicity)

plot(m_bty_gen$residuals ~ evals$gender)
abline(h = 0, lty = 3)  # adds a horizontal dashed line at y = 0

Exercise 8

The estimate of the beauty average changed to 0.07416 from 0.06664, but the R-squared is still low. So we can not be sure if it is due to chance.

Exercise 9

y_hat = bo_hat + b1_hat * bty_avg

Males tend to have the higher score.

multiLines(m_bty_gen)

Exercise 10

With three levels of rank, R shows three variable names and then (Intercept) for the fourth.

m_bty_rank <- lm(score ~ bty_avg - gender + rank, data = evals)
plot(m_bty_rank)

summary(m_bty_rank)

## 
## Call:
## lm(formula = score ~ bty_avg - gender + rank, data = evals)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.8713 -0.3642  0.1489  0.4103  0.9525 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       3.98155    0.09078  43.860  < 2e-16 ***
## bty_avg           0.06783    0.01655   4.098 4.92e-05 ***
## ranktenure track -0.16070    0.07395  -2.173   0.0303 *  
## ranktenured      -0.12623    0.06266  -2.014   0.0445 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5328 on 459 degrees of freedom
## Multiple R-squared:  0.04652,    Adjusted R-squared:  0.04029 
## F-statistic: 7.465 on 3 and 459 DF,  p-value: 6.88e-05

Exercise 11

I would expect the class level to not have much of an association with professor score since students of all levels feel different ways about professors.

m_full <- lm(score ~ rank + ethnicity + gender + language + age + cls_perc_eval 
             + cls_students + cls_level + cls_profs + cls_credits + bty_avg 
             + pic_outfit + pic_color, data = evals)
summary(m_full)

## 
## Call:
## lm(formula = score ~ rank + ethnicity + gender + language + age + 
##     cls_perc_eval + cls_students + cls_level + cls_profs + cls_credits + 
##     bty_avg + pic_outfit + pic_color, data = evals)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.77397 -0.32432  0.09067  0.35183  0.95036 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            4.0952141  0.2905277  14.096  < 2e-16 ***
## ranktenure track      -0.1475932  0.0820671  -1.798  0.07278 .  
## ranktenured           -0.0973378  0.0663296  -1.467  0.14295    
## ethnicitynot minority  0.1234929  0.0786273   1.571  0.11698    
## gendermale             0.2109481  0.0518230   4.071 5.54e-05 ***
## languagenon-english   -0.2298112  0.1113754  -2.063  0.03965 *  
## age                   -0.0090072  0.0031359  -2.872  0.00427 ** 
## cls_perc_eval          0.0053272  0.0015393   3.461  0.00059 ***
## cls_students           0.0004546  0.0003774   1.205  0.22896    
## cls_levelupper         0.0605140  0.0575617   1.051  0.29369    
## cls_profssingle       -0.0146619  0.0519885  -0.282  0.77806    
## cls_creditsone credit  0.5020432  0.1159388   4.330 1.84e-05 ***
## bty_avg                0.0400333  0.0175064   2.287  0.02267 *  
## pic_outfitnot formal  -0.1126817  0.0738800  -1.525  0.12792    
## pic_colorcolor        -0.2172630  0.0715021  -3.039  0.00252 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.498 on 448 degrees of freedom
## Multiple R-squared:  0.1871, Adjusted R-squared:  0.1617 
## F-statistic: 7.366 on 14 and 448 DF,  p-value: 6.552e-14

Exercise 12

The values such as R-squared and the summary p-value are very similar to the full model. Suggesting that class level really does not matter that much.

m_modified <- lm(score ~ rank + ethnicity + gender + language + age + cls_perc_eval 
             + cls_students + cls_profs + cls_credits + bty_avg 
             + pic_outfit + pic_color, data = evals)
summary(m_modified)

## 
## Call:
## lm(formula = score ~ rank + ethnicity + gender + language + age + 
##     cls_perc_eval + cls_students + cls_profs + cls_credits + 
##     bty_avg + pic_outfit + pic_color, data = evals)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.76617 -0.31818  0.09325  0.35169  0.93485 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            4.0938164  0.2905587  14.089  < 2e-16 ***
## ranktenure track      -0.1419977  0.0819039  -1.734 0.083656 .  
## ranktenured           -0.0895639  0.0659238  -1.359 0.174957    
## ethnicitynot minority  0.1383380  0.0773580   1.788 0.074405 .  
## gendermale             0.2046338  0.0514798   3.975  8.2e-05 ***
## languagenon-english   -0.2109231  0.1099296  -1.919 0.055655 .  
## age                   -0.0087375  0.0031258  -2.795 0.005407 ** 
## cls_perc_eval          0.0053940  0.0015382   3.507 0.000499 ***
## cls_students           0.0003430  0.0003622   0.947 0.344081    
## cls_profssingle       -0.0150776  0.0519931  -0.290 0.771956    
## cls_creditsone credit  0.4692494  0.1116766   4.202  3.2e-05 ***
## bty_avg                0.0412068  0.0174728   2.358 0.018785 *  
## pic_outfitnot formal  -0.1216791  0.0733912  -1.658 0.098026 .  
## pic_colorcolor        -0.1955247  0.0684550  -2.856 0.004486 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.498 on 449 degrees of freedom
## Multiple R-squared:  0.1851, Adjusted R-squared:  0.1615 
## F-statistic: 7.845 on 13 and 449 DF,  p-value: 3.699e-14

Exercise 13

The ethnicity changes raises score by 0.1234929 points, all other variables constant

Exercise 14.

The coefficients are mostly slightly smaller, with the signficance levels also being lower and higher.

m_modified2 <- lm(score ~ rank + ethnicity + gender + language + age + cls_perc_eval 
             + cls_students + cls_level + cls_credits + bty_avg 
             + pic_outfit + pic_color, data = evals)
summary(m_modified2)

## 
## Call:
## lm(formula = score ~ rank + ethnicity + gender + language + age + 
##     cls_perc_eval + cls_students + cls_level + cls_credits + 
##     bty_avg + pic_outfit + pic_color, data = evals)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.7836 -0.3257  0.0859  0.3513  0.9551 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            4.0872523  0.2888562  14.150  < 2e-16 ***
## ranktenure track      -0.1476746  0.0819824  -1.801 0.072327 .  
## ranktenured           -0.0973829  0.0662614  -1.470 0.142349    
## ethnicitynot minority  0.1274458  0.0772887   1.649 0.099856 .  
## gendermale             0.2101231  0.0516873   4.065 5.66e-05 ***
## languagenon-english   -0.2282894  0.1111305  -2.054 0.040530 *  
## age                   -0.0089992  0.0031326  -2.873 0.004262 ** 
## cls_perc_eval          0.0052888  0.0015317   3.453 0.000607 ***
## cls_students           0.0004687  0.0003737   1.254 0.210384    
## cls_levelupper         0.0606374  0.0575010   1.055 0.292200    
## cls_creditsone credit  0.5061196  0.1149163   4.404 1.33e-05 ***
## bty_avg                0.0398629  0.0174780   2.281 0.023032 *  
## pic_outfitnot formal  -0.1083227  0.0721711  -1.501 0.134080    
## pic_colorcolor        -0.2190527  0.0711469  -3.079 0.002205 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4974 on 449 degrees of freedom
## Multiple R-squared:  0.187,  Adjusted R-squared:  0.1634 
## F-statistic: 7.943 on 13 and 449 DF,  p-value: 2.336e-14

Exercise 15

y_hat = bo_hat + b1_hatethnicity + b2_hatgender + b3_hatlanguage + b4_hatage + b5_hatperc_eval + b6_hatcredits + b7_hatbty_avg + b8_hatpic_color

final_model <- lm(score ~ ethnicity + gender + language + age + cls_perc_eval + cls_credits + bty_avg + pic_color, data = evals)
summary(final_model)

## 
## Call:
## lm(formula = score ~ ethnicity + gender + language + age + cls_perc_eval + 
##     cls_credits + bty_avg + pic_color, data = evals)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.85320 -0.32394  0.09984  0.37930  0.93610 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            3.771922   0.232053  16.255  < 2e-16 ***
## ethnicitynot minority  0.167872   0.075275   2.230  0.02623 *  
## gendermale             0.207112   0.050135   4.131 4.30e-05 ***
## languagenon-english   -0.206178   0.103639  -1.989  0.04726 *  
## age                   -0.006046   0.002612  -2.315  0.02108 *  
## cls_perc_eval          0.004656   0.001435   3.244  0.00127 ** 
## cls_creditsone credit  0.505306   0.104119   4.853 1.67e-06 ***
## bty_avg                0.051069   0.016934   3.016  0.00271 ** 
## pic_colorcolor        -0.190579   0.067351  -2.830  0.00487 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4992 on 454 degrees of freedom
## Multiple R-squared:  0.1722, Adjusted R-squared:  0.1576 
## F-statistic:  11.8 on 8 and 454 DF,  p-value: 2.58e-15

Exercise 16

There are few to no outliers shown in the histogram. There are no obvious patterns in the residuals vs fitted data There are no obvious patterns viewing the data in order. Against both ethnicity and gender, the boxplots do not vary much.

hist(final_model$residuals)

plot(abs(final_model$residuals) ~ evals$score)

plot(final_model$residuals)

plot(final_model)

plot(final_model$residuals ~ evals$ethnicity)

plot(final_model$residuals ~ evals$gender)

Exercise 17

The fact that each of the observations came from only one university might make the samples not as independent as they should be for analysis using linear regression.

Exercsise 18

Male, non-minorities who teach a one-credit course, who had many students rate their course, received educations in English, receive high beauty scores, and are are young, are likely to have higher scores.

Exercise 19

I would not be comfortable generalizing these conclusions to professors at any university because these observations all came from one specific university in one small part of the world. It is likely that different cultural values around the world would make the results of this study difficult to generalize to all universities.